Data loading

AnalyzAIRR proposes two different data loading processes depending on the alignment files format.

l <- list.files(system.file(file.path('extdata/mixcr'),
                            package='RepSeq'),
                            full.names = T)

#> [1] "/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RepSeq/extdata/mixcr/tripod-30-813.txt.gz" 
#> [2] "/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RepSeq/extdata/mixcr/tripod-30-815.txt.gz" 
#> [3] "/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RepSeq/extdata/mixcr/tripod-31-846.txt.gz" 
#> [4] "/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RepSeq/extdata/mixcr/tripod-31-848.txt.gz" 
#> [5] "/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RepSeq/extdata/mixcr/tripod-35-970.txt.gz" 
#> [6] "/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RepSeq/extdata/mixcr/tripod-35-972.txt.gz" 
#> [7] "/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RepSeq/extdata/mixcr/tripod-36-1003.txt.gz"
#> [8] "/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RepSeq/extdata/mixcr/tripod-36-1005.txt.gz"

metadata <- read.table(system.file(file.path('extdata/sampledata.txt'),
                            package='RepSeq'), 
                            header=T,
                            row.names = 1)
metadata$cell_subset <- factor(metadata$cell_subset)
metadata$sex <- factor(metadata$sex)

RepSeqData <- readAIRRSet(fileList = l,
                          fileFormat = "MiXCR",
                          chain = "TRA",
                          sampleinfo = metadata,
                          keep.ambiguous = FALSE,
                          keep.unproductive=FALSE,
                          filter.singletons=FALSE,
                          aa.th=8,
                          outFiltered=TRUE,
                          raretab=FALSE,
                          cores=1L)
#> Running on 1 cores...
#> Loading and filtering sequences...
#> Assembling sequences...
#> Creating a RepSeqExperiment object...
#> Done.

The RepSeqExperiment object

assay(RepSeqData)

mData(RepSeqData)

oData(RepSeqData)$filtered

History(RepSeqData)

Filtering

• Filter sequences based on their count in each sample.

  filterCount() filters out sequences having a count below the chosen threshold. The function returns a RepSeqExperiment object that can be used to perform in-depth analyses on the remaining dataset

• Extract sequences that are either shared across samples, or private to a single sample.

  getPublic() allows to subset a RepSeqExperiment object in order to keep sequences that are shared by:

  • samples belonging to a specified group
  • samples within the whole dataset if the group parameter is not specified.
    
    

  The level of sharing can be specified with the prop parameter to:

  • keep sequences with a sharing level greater than or equal to the chosen threshold
  • keep sequences present in at least two samples if prop is not specified.
    

The function returns a RepSeqExperiment object with the selected clones.

# Get clones present in at least 50% of the samples belonging to the amTreg group
getPublic(x=RepSeqData, level="clone", group=c("cell_subset", "amTreg"), prop=0.5)

Similarly, it is possible to extract private sequences with getPrivate(). If the parameter singletons is set to TRUE, only private sequences with a count of 1 will be returned.

Note: All the above-mentioned functions can be applied at any of the following repertoire level: clone, clonotype, CDR3aa and CDR3nt.

• Filter productive or unproductive sequences

  getProductive() and getUnproductive() allow the extraction of productive and unproductive sequences respectively in case no filters were applied during the building of the RepSeqExperiment object.

• Exclude one or multiple samples from a RepSeqExperiment object.

  dropSample() offers the possibility to filter out one or multiple repertoires by specifying their corresponding sample_id. These repertoires can be ones identified as outliers during the exploratory analysis routine described further below.

Normalization

AnalyzAIRR proposes two normalization strategies:

• Down-sampling

This strategy can be applied when the studied samples largely differ in their repertoire sizes.

By using the sampleRepSeqExp() function, users can choose the value to which all the samples are downsampled. If not specified, the lowest number of sequences across the dataset will be used. This function returns a new RepSeqExperiment object with the downsized data.

RepSeqData_ds<- sampleRepSeqExp(x=RepSeqData, sample.size=50000)

• Shannon-based normalization

This strategy adapted in (Chaara et al., 2018) can be used to eliminate “uninformative” sequences resulting from experimental noise. It is particularly efficient when applied on small samples as it corrects altered count distributions caused by a high-sequencing depth.

The function ShannonNorm() allows the application of this strategy without the need to specify any parameter and returns a new RepSeqExperiment object with the corrected data.

RepSeqData_sh<- ShannonNorm(x=RepSeqData)

Eliminated sequences in each sample are stored in the otherData slot.

Basic statistics

• Visualize the metaData statistics

Statistics in the metaData slot, calculated during the building of the RepSeqExperiment object, can be visualized for each sample using the plotStatistics() function with the parameter groupBy set to NULL.

plotStatistics(x = RepSeqData, 
               stat = "nSequences",
               groupBy = NULL)

Figure 1: A barplot showing the number of TCR sequences identified in each sample.

plotStatistics(x = RepSeqData, 
               stat = "clonotype",
               groupBy = NULL)

Figure 2: A barplot showing the number of clonotypes identified in each sample.

• Calculate the count or the proportion of a chosen repertoire level.

  countFeatures() returns the calculated values for all the samples within the RepSeqExperiment object.

  The function takes into account the weight of the studied level, i.e. the number of sequences expressing a specific gene segment, or the count of a sequence in a sample.

countFeatures(x=RepSeqData, level="J", scale="frequency")

Diversity estimation

• Rarefaction analysis

  Species richness can be assessed by a rarefaction table generated using the rarefactionTab() function.

  Rarefaction curves can as well be plotted using plotRarefaction()

plotRarefaction(x=RepSeqData, colorBy="sample_id")
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.

Figure 3: Rarefaction curves plotting the number of clonotypes against the number of sequences in a sample.

• Calculation of diversity indices.

  The function diversityIndices() computes a set of diversity indices on a chosen repertoire level for each sample. The calculated indices are the following:

  • Shannon index: Calculates the proportional abundance of species in a repertoire (Shannon, 1948).

  • Simpson index: a measure taking into account the number of species present as well as their abundance. It gives relatively little weight to the rare species and more weight to the frequent ones (Simpson, 1949).

  • Inverse Simpson index: Is the effective number of species that is obtained when the weighted arithmetic mean is used to quantify average proportional abundance of species.

  • Gini coefficient: Measures the degree of inequality in a distribution of abundances (Gini, 1921).

  • Chao1: Estimates undetected species using the information on the rarest species (the numbers of singletons and doubletons) (Chao, 1984).

  • Improved Chao1: An extension of Chao1 which uses additional information, namely, the numbers of tripletons and quadrupletons (Chiu et al., 2014).

diversityIndices(x=RepSeqData, level="clone")

The calculated indices can be plotted using plotDiversity(). The parameter index allows to choose which index is to be plotted.

If groupBy is not specified, calculated values in each sample are represented.

plotDiversity(x=RepSeqData, 
              level = "clone", 
              groupBy = NULL, 
              index = "gini")

Figure 4: A barplot showing the Gini values calculated at the clone level within each sample.

plotDiversity(x=RepSeqData, 
              level = "clonotype", 
              groupBy = NULL, 
              index = "shannon")

Figure 5: A barplot showing the Shannon values calculated at the clonotype level within each sample.

• Evaluation of the Renyi diversity

  The Renyi index is a generalization of the Shannon entropy. It represents the distribution of clonal expansions as a function of the parameter alpha. It equally considers all species including the rare ones at alpha=0, and up-weighs the abundant species with an increasing value of alpha. Alpha =1 is an approximation of the Shannon index; alpha = 2 corresponds to the Simpson index and alpha=Inf corresponds to the Berger-Parker index. The latter highlights the highest clonal expansion in a repertoire.

  Renyi values can be computed at any repertoire level using renyiIndex(). By default, it is assessed for the following alpha values: \[\alpha \in {0, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, Inf}\]

renyiIndex(x=RepSeqData, level="clone")

Similarly to the previous function, Renyi’s entropy values can be plotted using plotRenyiIndex().

plotRenyiIndex(x=RepSeqData, level="J", grouped=FALSE, colorBy = "sample_id")

Figure 6: The Renyi values were calculated at the J gene level for each of the pre-defined alpha values and plotted for each sample.

Clonal distribution

• Clonal distribution within count intervals

  The plotCountIntervals() function can be used to assess the clonal distribution within pre-defined count intervals.

  It computes:

  • The proportion that a defined repertoire level within each count interval represents in the total repertoire, thus giving a global evaluation of the count distribution.
  • The cumulative frequency of a defined repertoire level within each count interval thus giving the proportion that each interval occupies in the repertoire.

plotCountIntervals(x=RepSeqData, level="clone")

Figure 7: Stacked bars showing the proportion (left plot) and the cumulative frequency (right plot) of clones within each count interval in each sample.

• Clonal Distribution per decreasing rank

Another way to assess the clonal distribution is by visualizing the distribution as a function of the ranks. Samples can be colored by sample_id as well as by any group column from the mData() slot.

plotRankDistrib(x=RepSeqData,level="clone", colorBy="sample_id", grouped=F, scale="frequency")

Figure 8: Plotted is the distribution of the clones as a function of their occurrence rank within each sample.

Basic statistics comparison

The majority of the functions described in this section can be applied to compare groups of samples and perform statistical tests.

• Compare metaData statistics

  The calculated statistics in the metaData slot can be compared between groups of samples by specifying the groupBy parameter in the plotStatistics() function.

plotStatistics(x = RepSeqData, 
               stat = "clonotype", 
               groupBy = "cell_subset")

Figure 15: The number of sequences is not statistically different between the two studied cell subsets. Boxplots represent the median across all samples belonging to the same group. A Wilcoxon test is applied and adjusted p-values using the Holm method are shown.

plotStatistics(x = RepSeqData, 
               stat = "clonotype", 
               groupBy = "sex")

Figure 16: The number of sequences is statistically different between the two studied cell subsets. Boxplots represent the median across all samples belonging to the same group. A Wilcoxon test is applied and adjusted p-values using the Holm method are shown.

• Compare repertoire diversities

  The plotDiversity() function also supports the groupBy parameter to compare groups for a chosen diversity index.

  For the Renyi index, the grouped parameter should be set to TRUE.

plotDiversity(x = RepSeqData, 
              level = "clone", 
              groupBy = "cell_subset", 
              index = "gini")

Figure 17: The gini index is calculated at the clone level and compared between cell subsets. Gini values are statistically different between amTregs and nTregs. Boxplots represent the median across all samples belonging to the same group. A Wilcoxon test is applied and adjusted p-values using the Holm method are shown.

plotRenyiIndex(x = RepSeqData, 
                  level = "clone", 
                  grouped = "TRUE", 
                  colorBy = "cell_subset")

Figure 18: The Renyi values were calculated at the clone level for each of the pre-defined alpha values and plotted for each sample. Groups represented in differing colors. For each group, the mean (circles) and stantard error (shade) are plotted.

• Compare clonal distributions

  Similarly to the analysis applied at the single-sample level, plotCountIntervals() allows to plot and compare clonal distributions within defined count intervals between experimental groups.

plotCountIntervals(x=RepSeqData, level="clone", groupBy = "cell_subset")

Figure 19: The proportion (left plot) and the cumulative frequency (right plot) of clones within each count interval was compared between amTreg and nTreg samples. A Wilcoxon test is applied and adjusted p-values using the Holm method are shown.

Similarity analysis

The repertoire sharing at any level evaluates the degree of convergence between repertoires and experimental conditions. Different statistical methods are proposed herein to evaluate this convergence.

• Compute the number of shared sequences

  The number of shared sequences, at any level, between samples belonging for instance to the same experimental group can be obtained using the plotEulerr() function. If sampleNames is not specified, the first 3 samples in the datasets will be analyzed.

ctrnames <- rownames(mData(RepSeqData))[which(mData(RepSeqData)[,"cell_subset" ] %in% "nTreg")]
plotEulerr(x = RepSeqData, 
           level = "clone",
           sampleNames = ctrnames)

Figure 20: Euler diagram representing the number of shared clones between the nTreg samples.

• Assess the correlation between a pair of samples

  The correlation between a pair of repertoires can be calculated using the plotScatter() function by simply specifying a pair of sampleNames to compare.

plotScatter(x = RepSeqData, 
            level = "V",
            scale = "frequency", 
            sampleNames = c("tripod-30-813","tripod-31-846"))

Figure 21: A biplot comparing the V gene usage between two defined samples. A linear model is fitted on the data. The V usage is significantly similar between the two samples.

plotScatter(x = RepSeqData, 
            level = "clone", 
            scale = "frequency", 
            sampleNames = c("tripod-30-813","tripod-31-846"))

Figure 22: A biplot comparing the clonal expression between two defined samples. A linear model is fitted on the data. Clonal overlap and expression are significantly low between the compared samples.

• Compute dissimilarity distances

  RepSeq proposes a list of dissimilarity indices, each taking into account different parameters. The proposed methods include:

  • The Jaccard similarity: a measure of similarity between sample sets defined as the size of the intersection divided by the size of the union of the sample sets.

  • The Morisita-Horn similarity: a measure of similarity that tends to be over-sensitive to abundant species.

  Details on these indices and others can be found here

  These distances can be calculated at any level of the repertoire and can be:

  • Visualized as a dissimilarity heatmap using plotDissimilarityMatrix(). This function performs a hierarchical clustering on the calculated distance scores using the method specified in the clustering parameter.
  • Used to perform a multidimensional scaling analysis or a principal component analysis using plotDimReduction() with the dim_method parameter set to MDS

plotDissimilarityMatrix(x = RepSeqData, 
                        level = "clone", 
                        method = "morisita", 
                        clustering = "ward.D",
                        binary = "FALSE")

Figure 23: A heatmap showing the Morisita distance calculated at the clone level between all samples within the dataset. The ward.D method was used to perform a hierarchical clustering on the calculated distances.

plotDimReduction(x = RepSeqData, 
                 level = "clone",
                 method = "morisita", 
                 colorBy = "cell_subset", 
                 dim_method = "MDS")

Figure 24: The Morisita distance were calculated at the clone level between all samples within the dataset and used to perform a multidimensional scaling analysis.

Differential analysis

Differentially expressed genes or sequences can be identified using diffExpGroup. The experimental groups to be compared can be specified with the group parameter, and the function outputs the statistics calculated for each gene/sequence.

In the following example, we identify over-expressed VJ combinations within the amTreg samples compared to the nTregs.

DS <- diffExpGroup(x = RepSeqData,
                   colGrp = "cell_subset" ,
                   level = "VJ",
                   group = c("cell_subset", "amTreg", "nTreg"))